Even-odd effects in NSN scattering problems: Application to graphene nanoribbons 
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We study crossed Andreev reflection (CAR) of electrons or holes in normal metal-superconductor- 
normal metal junctions and highlight some very strong effects of the underlying lattice. In particular, 
we demonstrate that for sharp interfaces and under certain, albeit generic, symmetry conditions, the 
CAR probability exactly vanishes for an even number of atoms in the superconducting region. This 
even-odd effect applies notably to NSN junctions made of graphene nano-ribbons with armchair 
edges and for zigzag edges with somewhat more restrictive conditions. We analyze its robustness 
towards smoothing of the boundaries or doping of the sample. 
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I. INTRODUCTION 

A significant amount of work has recently been devoted 
to the study of possible realizations of spin-entanglers 
in solid-state devices. Andreev reflection 1 at the inter- 
face between a normal metal and a superconductor is 
of particular interest in this context. Indeed the con- 
version of an incoming hole into a reflected electron at 
the interface is readily identified with the injection of a 
Cooper pair into the normal lead and its splitting into 
two spin-entangled electrons. 2 When two normal leads 
are connected to a superconducting barrier, a non-local 
electron-hole conversion 3 ' 4 can take place, as a Cooper 
pair is splitted among both leads. Such a non-local pro- 
cess could be enhanced by Coulomb interaction 5-11 or by 
energy-filtering processes 12,13 . Non-local or crossed An- 
dreev reflection (CAR) in normal metal-superconductor- 
normal metal (NSN) junctions was observed experimen- 
tally, 14,15 and shown to depend strongly on the thickness 
of the superconductor. CAR is in general not perfect 
and will coexist with cotunneling, that is, the transfer of 
an electron or a hole across the superconducting barrier, 
as well as local normal or Andreev reflection. In non- 
local conductance measurements, one measures the sum 
of both effects and it is therefore not possible to measure 
the CAR signal alone. However, more recently, spatial 
separation of the CAR signal from the electron cotunnel- 
ing has been proposed theoretically 16,17 . In setups where 
a bias voltage only exists between the superconductor 
and the normal metal leads, the electron cotunneling does 
not appear and one can use current noise measurements 
to distinguish CAR from the local Andreev reflection 18,19 
as well as conductance measurements in more involved se- 
tups containing quantum dots 20,21 in the Coulomb block- 
ade regime between the superconductor and the normal 
metal leads or Aharonov-Bohm rings 10 ' 22 . 

NSN junctions could in principle be made of graphene 
nanoribbons, 12 ' 2 where superconductivity is induced 
by the proximity to a nearby superconductor. 26-29 At 
interfaces between graphene and superconductors, new 
physics emerges because of the peculiar band structure 
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FIG. 1: NSN junction made from a graphene nanoribbon 
with armchair edges. In the S region, an onsite pairing poten- 
tial A is induced by proximity to an s-wave superconductor. 

of the carbon monolayer. For instance, in a single NS 
junction, specular local Andreev reflection was predicted 
to exist, in striking contrast with conventional metal 
leads, 30 which has observable effects in the transport 
characteristics 22 ' 30-32 . In graphene NSN junctions, it 
appears to be possible to favor CAR over co-tunneling 
by either considering asymmetrically doped 12 or fer- 
romagnetic leads. 23,25 Interestingly, as the length of 
the superconductor is increased, the CAR probability 
exhibits oscillations, a consequence of quantum in- 
terferences in the superconducting region 25 , a feature 
shared by all NSN junctions with interfaces sharp on the 
scale of the Fermi-wave length in the superconductor. 
In our paper, we carefully study these oscillations 
and go beyond by looking at length dependence on 
the scale of the lattice constant. In particular, we dis- 
cover an effect related to the parity of the S region length. 

This effect is in fact not graphene-specific but appears 
already in the simplest model of an NSN junction, namely 
a one-dimensional lattice with nearest-neighbor hopping, 
where an on-site pairing potential varies from zero in the 
normal leads to a value A ^ in the superconducting re- 
gion. If the NS interfaces are atomically sharp, as shown 
in Fig. 3 below, and if the following two conditions are 
fulfilled 

(i) the spectrum has particle-hole symmetry, that is, 
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the Fermi energy is zero everywhere in the sample 
and 

(ii) the superconducting region spans an even number 
of sites, 

then the CAR probability vanishes exactly. For an odd 
number of sites in the S region, the CAR probability does 
not vanish. To the best of our knowledge, such a strong 
even-odd effect 42 was not reported before. This effect 
carries over to armchair graphene nanoribbons, if con- 
dition (ii) is interpreted as requiring an even number of 
transverse zigzag chains in the S region. This situation 
is represented in Fig. 1 and the CAR probability oscil- 
lations, as obtained from wave-function matching, are 
plotted in Fig. 2. In armchair nano-ribbons, transverse 
modes are completely decoupled and their contribution to 
the CAR signal add up incoherently. It turns out that all 
transverse modes exhibit the even-odd effect under con- 
ditions (i) and (ii) . Therefore, the CAR probability actu- 
ally vanishes for even lengths of the S region, independent 
of the excitation energy, below the gap A. The zeroes of 
the CAR probability are analytically understood from an 
analysis of the symmetries of the problem, at the level of 
the transfer matrix. In zigzag ribbons the situation is 
complicated by the strong dependence of Andreev reflec- 
tion on the ribbon width. 33 Interestingly, we found the 
same even-odd effect for anti-zigzag ribbons, yet only for 
the lowest mode. Strikingly, under the symmetry condi- 
tions discussed above, the lattice has a strong influence 
on the transport characteristics, even when the super- 
conducting gap evolves smoothly, but symmetrically, at 
the two NS-interfaces. Note that such a behavior is not 
observed in the Bogoliubov-de Gennes equation based 
on the Dirac equation 30 where CAR should vanish when 
the Fermi energy is zero due to pseudospin conservation. 
We observe indeed a very small CAR signal in the tight- 
binding model, however, the even odd effect remains im- 
portant. The total signal of the CAR could be enhanced 
significantly in a setup with many modes (note the ef- 
fect is observed in any armchair mode), this is realized 
in particular when the ribbon width is large. In the case 
of the zigzag ribbon, the transverse modes will couple in 
a scattering experiment and the situation is more subtle. 

The article is organized as follows. In Section II, we 
first demonstrate the even-odd effect for the ID lattice. 
In particular, we show that the even-odd effect is in- 
dependent of the excitation energy, or the strength of 
the pairing potential. Subsequently, an extension to 2D 
square lattices is presented. Then, in Section III, we 
transpose the analysis to graphene nano-ribbons. Finally, 
in Section IV, we discuss non-idealities, in particular dop- 
ing in the S region and smoothing of the NS interfaces, 
and interpret numerical results for the CAR probability, 
obtained with the recursive Green's function method, in 
the light of the even-odd effect. Numerical results for the 
zigzag-ribbon are analysed in the last subsection. De- 
tailed calculations can be found in the appendices. 
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FIG. 2: (Color online) Crossed Andreev reflection (CAR) 
probability for the NSN junction of Fig 1, with W = 92, 
the number of dimers in the transverse direction; the ribbon 
is metallic. The Fermi level is at the Dirac point in the whole 
sample. The excitation energy is e = 0.03t, and the pairing 
gap A = O.lt, with t the nearest-neighbor hopping amplitude. 
Only the metallic mode is excited. The CAR signal vanishes 
exactly for even N. 



II. EVEN-ODD EFFECT IN THE LINEAR 
CHAIN 
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FIG. 3: NSN junction in a ID infinite lattice. The S region 
(shaded box) extends over L sites. 

We first consider a NSN junction in the simplest one- 
band model: a chain with nearest-neighbor hopping. The 
position of a given site is labeled by an integer i running 
from — co to +oo and the S region extends over L sites. 
We choose, without loss of generality, sites I — \ and 
I = L to be the end-sites of the S region, as exemplified 
in Fig. 3. Electrons moving on such a lattice, with lat- 
tice constant a, are then described by the tight-binding 
Hamiltonian 



H 



h 



H.c 



(1) 



where m,n run over all integers, a =t,i is a spin 1/2 
index and h mn = —t5 njm ±i describes nearest-neighbor 
hopping while A m?l = A5 m _ n , for 1 < m < L, is the 
pairing potential in the S region. Solving the scattering 
problem for a given excitation energy e, that is, find- 
ing the scattering matrix S(e), amounts to solving the 
Bogoliubov-de Gennes (BdG) equations 



= e 



(2) 
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where 



'HBdG = 



n -E F 

A 



A 

-Ho + E F 



(3) 



with [Ho]mn = Knn and [A]„ m = A mn . E F is the Fermi 
energy. We have dropped the spin index due to spin- 
rotation invariance of the Hamiltonian and chosen A to 
be real, (resp. ^h) describes electron-like (resp. hole- 
like) excitations. We will look for solutions of the trans- 
lational invariant problem corresponding to each region 
(N or S) and match the wave-functions at the NS bound- 
aries. 



A. Scattering states and transfer matrix 




FIG. 4: (Color online) Spectrum for electron-like (blue, solid) 
and hole-like (red, dashed) excitations, for Ef = 0. We have 
set t = a = 1. For a given excitation energy e (dotted hori- 
zontal line) there are four allowed solutions which correspond 
to right and left-moving electrons or holes. 

In the normal regions, where A mn = 0, a general so- 
lution of Eq. (2), for a given excitation energy e, is a 
superposition of incoming and outgoing plane waves, 
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where we have defined, for electron-like excitations, right- 
going waves (+) and left-going waves (— ) as "$>f e = 
e ±ik c £a wn ji e f or hole-like excitations, right-going waves 
and left-going waves are of the form ~ e =Rfch£a . 

The wave-vectors k e and kh are the positive roots of 
e = £,(k) and e = — respectively, with £(fc) = 
— 2icos(fcet) — Ep (see Fig. 4). Particle-hole symme- 
try of the spectrum, for Ep = 0, implies the important 
constraint 



n/a, 



(5) 



for any excitation energy e. The amplitudes in Eq. (4) 



(aj,e,aj,h) , Bj 



for j = 1,2 



are related by the scattering matrix S through 

= S 



Br 
B 2 



Ai 
A 2 



(6) 



The scattering matrix is here a 4 x 4 matrix which has 
the generic form 



r t' 
t r' 



(7) 



with r, t, t', r', 2x2 matrices. We find it more convenient 
to work with transfer matrices - so-called M-matrices - 
which relate amplitudes on the right of the scattering 
potential to amplitudes on the left, that is 



B 2 
A 2 



We also introduce a block notation for M as 



M = 



7 



(8) 



(9) 



where, again, a, (3, 7 and S are 2x2 matrices. Some im- 
portant properties of transfer matrices are summarized in 
Appendix A. In particular, current conservation implies 
that 



t = 



(10) 



We show in the next sections that, under the two condi- 
tions stated in the introduction, the off-diagonal elements 
of t are zero, and therefore the CAR probability vanishes. 



B. Wave- matching at the NS boundaries 

In practice, the transfer matrix as well as the scattering 
matrix are obtained by searching for plane- wave (possibly 
with a complex wave-vector) solutions in the supercon- 
ducting region and matching the wave-functions at the 
two boundaries. In the S region, plane waves ^s,£ are 
solutions of 



A(k) A 
A -A(fc) 



(11) 



with A(fc) = —2tcos(ka) — Ep — Vs, and Vs the doping in 
the S region. Eigenvalues are e — ±y / A(k) 2 + A 2 , with 
eigenvectors X±W = (F±(k), 1) T and 



F ± (k) - (A(fc) ± ^A(fc) 2 + A 2 )/A. (12) 

There are, generically, four solutions for a given excita- 
tion energy e. For |e| > A, wave-vectors are real, with 
two right and two left-going waves. However, we are in- 
terested in sub-gap transport, that is, \e\ < A, for which 
electron excitations can be converted into hole excita- 
tions, and vice-versa, via the creation or the annihilation 
of a Cooper pair. In this case, the solutions are four 
evanescent waves with complex wave-vectors k given by 



,1 E F + V S ± WA 2 - e 2 

k = ±-arccos 

a \ 2t 



(13) 
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Two of these wave- vectors have a positive imaginary part 
and correspond to right-decaying waves while the other 
two have a negative imaginary part and correspond to 
left-decaying waves. We define kg as the wave-vector 
with positive real and imaginary parts. From Eq. (13) it 
is clear that — k* s is the other right-decaying mode while 
the left-decaying modes are given by — k$ and kg. A 
general solution in the S region is of the form 



(a) 



+ a s .iXa{-k s )e~ lksta + a s ,2Xa{k* s )e ik ^ a , (14) 



with <j = sgn(e). For a particle-hole symmetric spectrum 
(Ep + Vs = 0) we find the following constraint 
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FIG. 5: Decomposition of the NSN junction of Fig. 3 into 2 
successive NS and SN step junctions represented in diagrams 

(a) and (c) respectively. Note that in the S region, for |e| < A, 
the arrows indicate right or left decaying modes. Diagrams 

(b) and (c) are related by a translation of L lattice sites. 



<Rc(k s a) = -. (15) 
The wave-matching conditions 34 read 

(fc) = •« - fe) = 1„ (i6) 

at the left boundary between sites and 1, as well as 

fc)=*"» d fc:0=**- <i7 > 

at the right boundary between sites L and L + l. A 
(numerical) resolution of this system of 8 equations will 
yield all elements of the transfer or scattering matrices. 
In the next section, using a suitable decomposition of 
the transfer matrix, we show analytically, and without 
solving explicitly the system of equations, how the even- 
odd effect in the CAR probability is a direct consequence 
of the symmetries of the problem. 



C. Decomposition of the transfer matrix 

One advantage of working with transfer matrices 
is that upon decomposing the scattering region into 
several, successive, scattering elements, the total transfer 
matrix is nothing but the product of transfer matrices 
for each scattering element. For instance, the NSN 
junction of Fig. 3 is decomposed into two NS or SN 
steps, as represented in Fig. 5. 



The transfer matrix M of the NSN junction is the prod- 
uct of two transfer matrices, one for each step 

M = Mi(L)Mj;(0). (18) 

In Fig. 5 we have recasted the amplitudes introduced 
in Eq. (14) into A s = (a s ,i, a sa ) T , &s = Q>s,l, bsaV ■ 
Thus, the transfer matrices M#(0) and Mi(I) relate 
amplitudes in the S region to amplitudes on the left and 
on the right, respectively, through the equations 

(X)= m *<»>(b:). < is » 

and 

(SH^ft)- (2o) 

Equation (18) is a direct consequence of Eqs. (8), (19) 
and (20). Note that Eqs. (19), (20) are only reformula- 
tions of the wave-matching conditions Eqs. (16), (17), 
respectively. Interestingly, there is a simple relation be- 
tween Mjj(O) and Mi(L). First note that the junctions 
(c) and (b) in Fig. 5 are related by a translation of L lat- 
tice sites. Therefore, amplitudes for the two associated 
scattering problems only differ by phase factors, that is 

M L (L) = Dia.g[e- ik ° La ,e ikhLa ,e ik ° La ,e- ikhLa ] 
M i (0)Diag[e afesiQ , e~ ik 's La , e~ iksLa , e lk "s La }. (21) 

Next we work out the relation between M^(0) and 
Mfl(0). Due to the left-right symmetry of the problem, 
we also have 

M L (0) = [M^(O)]- 1 . (22) 

As we did for M, we introduce a block decomposition for 
M L (0) and M R (0) through 

MW0)= ( aR ^ (23) 

1 n ' \1R(L) $R(L)J V ^ 
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Following Eq. (10), we compute the upper-left block a From Eqs. (18), (22), (23) and using blockwise inver- 

in M. We factor out the free propagation and introduce sion formulae, a. L = [a R - PrS' 1 ^]' 1 and (3 L = 

-ik e La q \ - [cxr - PrS^jr] PrSr 1 , we arrive at 
e ik h La)oc- (24) 



J2k s La 





-i2SRe(fc s )La 



t R 1 PR S R 1 f o e i2SRc(fe s )ia) lR 



[1 2 - a^^R^Vfl] 



(25) 



So far, we have only used the left-right symmetry of the 
barrier, and Eq. (25) is very general. To prove our state- 
ment about the CAR probability we now simplify the 
expression for a -1 using the conditions (i) Ep = Vs = 
and (ii) L even, stated earlier. Recalling that particle- 
hole symmetry of the spectrum implies *Ke(fcsa) = ir/2 
— see Eq. (15) — it is now appropriate to use the de- 
composition 



1 

e i7rL 



1 



l 



AttL 



(26) 



from which it follows that if L is an even integer then 

& -l = e ik s La ^ e i2k s La _ + j^j - 1 > ^7) 

where we have introduced the matrix 

A=l 2 - {ol-rPr) (^Sfl) = I2 r' L r R . (28) 

In the last equation, tr and r' L are two unitary matrices, 
that can be identified with the reflection matrices of the 
NS steps of Fig 5(a) and Fig 5(b), respectively, see App. 
A. To prove that the CAR amplitude is zero, it is enough 
to prove that the matrix A has zero off-diagonal elements. 
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Mj(0) A ^ 
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Mi(0) 



FIG. 6: The NS junction in (b') is the same as the one in (b), 
only shifted by one lattice site to the left. Is is related to the 
junction in (a) by a parity transformation. 

The two matrices tr and r' L are related by a trans- 
lation of one lattice site. This is illustrated on Fig. 6, 
where we represented the NS junction (b') which is re- 
lated to (b) by a translation of one site and to (a) by a 
parity transformation - see App. A for a detailed proof. 



As a consequence, 



A = l 2 



e 



,—ik e a q 



Za)r R { e Q / kha ) TR. (29) 



Particle-hole symmetry of the spectrum also implies r# = 
Vx v r&x- It follows that off-diagonal elements of A are 
zero. Indeed we have: 

A eh _ e -2rk a a A + e r(k h + k c )a^ ^ee = ( 3 q) 

and 

jhe = e 2ik h a /j + e -t{k h + k<! )a^ ^hh = ^ ^ 

due to the constraint in Eq. (5). 

D. Extension to the 2D square lattice 

The extension of the present study to the 2D square 
lattice is straightforward. We consider an NSN junction 
on a ribbon of width W. a is again the lattice constant 
and t is the hopping energy in both longitudinal and 
transverse directions. Imposing open boundary condi- 
tions, the transverse wave- vector p is quantized as 



P 



(W + l)a 



with r = 1, 2, 3, . 



(32) 



and is independent of the longitudinal wave-vector k. 
The 2D scattering problem decomposes into several in- 
dependent ID scattering problems, one for each available 
transverse mode at a given excitation energy e. The spec- 
trum of the ribbon is 

£(fc,p) = -E F - 2tcos(pa) - 2tcos(ka). (33) 

This means that each transverse mode has an effective 
Fermi energy E' F — Ep + 2tcos(pa), that can in principle 
be tuned to zero. For such modes, with E' F = 0, one 
will observe the same even-odd effect as in the linear 
chain, namely the absence of CAR for even lengths of the 
superconducting region. But for an arbitrary excitation 
energy, many modes with E' F =^ are likely to be excited 
- this number grows with the width of the ribbon - and 
the even-odd effect will be obscured. 
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In the next section, we will see that armchair graphcne 
nanoribbons, due to their peculiar (sub-)lattice geometry, 
exhibit an even-odd effect at zero Fermi energy Ep, for 
each and every transverse mode. 



III. EVEN-ODD EFFECT IN GRAPHENE 
NANORIBBONS 

A. Armchair edges 

1. Bogoliubov-de Gennes equations 

We now consider a graphene nanoribbon with arm- 
chair edges, as represented in Figs. 7 and 8. Let t be 
the hopping energy between neighboring lattice sites, 
a the distance between sites on a given sublattice and 
ot = v3a the length of a unit-cell. The width of the 
ribbon is W and counts the number of dimer lines in 
the transverse direction. We follow closely the notation 
of Ref. 35. Let us write down the tight binding Hamil- 
tonian H = Hq^ + Hq^ + iJ pa i r for armchair graphene 
ribbons, with 



m odd 



* a l,m,a^e.m+l, 
i m 

+ H.c. 



^e,m,a a £+l,m+l,cr 



(34) 



for a =t, I and 



pair 



A 



I'm 



H.c. 



Here, a 



and b 



t 



(35) 

are creation operators on the A 

and B sublattices, respectively, while a lma and b lmcr 
are the conjugate annihilation operators. We look for 
plane-wave solutions in the normal region, by taking the 
Fourier transform of the Hamiltonian in the y direction. 
Furthermore, spin-rotation invariance of the problem al- 
lows us to drop the spin 1/2 index. Plane- wave solutions 
can thus be written as 



where (ipA(k), ipB(k)) T is a spinor of the form 



(p A (k) 



_ se -iO(k) 

1 



(36) 



(37) 



and the phase 8(k) is defined in Eq. (42) below. Several 
comments are in order. Due to open boundary conditions 
on the armchair edges, the transverse momentum p is 
quantized as 



2rir 



P 



(W + l)a 



with r = 1, 2, 3, . 



(38) 



We furthermore use the convention that 

Dl,\A = Vi.2B = Vi,3A = ■■■ = £a T , (39) 
Ve.iB = Ve,2A = Vi,3B = ■■■ = (£+ l/2)a T , (40) 

as indicated in Fig. 7. In other words we absorb the 
differences of coordinates between sites mA and (m+l)B 
in the phase 0(k). Spinors in Eq. (37) are eigenstates of 
the effective 2-band Hamiltonian 



e + e ~ ikaT / 2 



(41) 



with eigenvalues E s (k : p) = s\e p + e — lka T/2^ \Y e have 
defined e p = 2 cos(pa/2) and s = ±1 is a band index. One 
can check that the phase difference between the spinor 
components is given by 



, f i p -ika T /2 \ V2 



(42) 



Note that, including the pairing potential, the 
Bogolubov-de Gennes equations reduce to the following 
effective form 



n (k,p)-E F Al 2 \ (9 e 

Ala -n (k )P ) + E F J Uf h 



(43) 

and can be solved independently for each value of the 
transverse momentum p, at a given excitation energy e. 
^ e and are now spinors, and Ho(k) is the effective 
2-band Hamiltonian introduced in Eq. (41). 



2. Spectrum and scattering states 

The superconducting region where A is non-zero ex- 
tends over L unit-cells, in the sense that the left NS in- 
terface lies in the cell I = while the right interface lies 
in the cell £ = L, as pictured in Fig. 8. Anticipating what 
follows and building on the formalism introduced for the 
treatment of the linear chain, we write the transfer ma- 
trix M( iJ ) of such an NSN junction as a product of two 
simpler transfer matrices, M^'(O) and M^p(£), one for 
each NS interface. Thus, 



M (i,j) =M ( g ) (L)M ( £(Q), 



(44) 



The labels i and j can, in principle, take the values 
1,2,3,4 corresponding to four distinct positions of the 
interfaces, as depicted in Fig. 8. However only positions 

1 and 2 are inequivalent, as 3 and 1 are related by a 
translation of half the unit-cell; the same is true for 

2 and 4. Due to the two sublattices in graphene, 1 
and 2 are not related by a simple translation in the 
longitudinal direction. In this section, we restrict i and 
j to the values 1 and 2. There are in total, for fixed 
L 7 four a priori inequivalent configurations of the NSN 
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m — 1 

771= 2 



P 




• /I sublattice 
O B sublattice 



m = W 



FIG. 7: Unit-cell of an armchair graphene nanoribbon of 
width W. An atom has coordinates (ma/2,y = ye, m A(B)), 
where I is the index of the unit-cell and A, B is the sublattice 
index. 



junction, which are labeled by (i, j) — (1, 1), (2, 2), (1, 2), 
or (2,1). The condition of an even number of sites 
in the S region here translates to = (1,1) or 

(2,2), which actually corresponds to an even number 
of zigzag columns. We will show that under conditions 
(i) particle-hole symmetry of the spectrum, that is, 
Ep = Vs = and (ii) an even number of sites in the S 
region, the CAR probability vanishes, for each transverse 
mode. Incidently, if condition (i) is fulfilled then the 
odd configurations (1,2) and (2,1) are equivalent, a 
consequence of time-reversal symmetry. For these two 
configurations, the CAR probability is finite. 

First, let us define the scattering states and state the 
wave-matching conditions. We work for a given trans- 
verse mode, with momentum p. Similar to the ID lat- 
tice, we introduce scattering states as superpositions of 
incoming and outgoing plane-waves according to 



in the left region, and 



(45) 



*h,t,m = KhK,Lm + «2,^ M , m , (46) 

in the right region. The right (+) and left (— ) going 
wave-functions are defined in App. B. Their wave-vectors 
are given by the solutions of e — —Ep + E s (k, p) and e = 
Ep — E s '(k,p), for electron-like and hole-like excitations, 




' = 



FIG. 8: NSN junction on an armchair graphene nanoribbon. 
The S region spans L unit-cells, in the sense that the left 
boundary is in the th unit-cell while the right boundary is 
in the L th unit-cell. In each of these two cells there are two 
inequivalent positions for the NS junction, 1 and 2 or 3 and 
4 (see main text). 



respectively. We define k e and kh as the positive roots of 
these equations. An example is given in Fig. 9. Particle- 
hole symmetry of the spectrum, Ep = 0, implies 



kp — ku 



(47) 
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FIG. 9: (Color online) Electron-like (blue) and hole-like (red) 
excitations for the metallic mode of armchair graphene (W = 
32), with E F = 0.2. For s = 0.4 (horizontal dotted line), 
electron excitation are in the conduction band (solid line, s = 
1) and holes in the valence band (dashed line, s' = —1). 

In the S region, we take into account possible doping 
that shifts the Fermi energy by Vs. The BdG Hamil- 
tonian has there four eigenvalues and eigenvectors. For 
a given excitation energy, there are again four possible 
solutions. As far as sub-gap transport is concerned, we 
write the corresponding wave- vectors as k$, —k* s (right- 
decaying modes) and — ks, k$ (left-decaying modes) 
where ks is the solution with positive real and imaginary 
parts (see App. B for more details). In the particle-hole 
symmetric case, Ep + Vs = 0, the wave-vectors satisfy 



the important constraint 
^e(fcs) = 



if e p < 0, 
2ir/a.T if e p > 0. 



(48) 



A general solution of the BdG equations in the S region 
is a superposition of right and left decaying modes of the 
form 

*SAm = b s ,iT){ks)xi+(ks) + b s ,2®(-k* s )xi+(-ks) 

+ a s ,i®(-k s )xi+(-ks) + a s ,2®(k* s )xi+(k*s), 

(49) 

where we have introduced the diagonal matrix D(k) = 



Diag[e 



iky t 



a iky t 



and Xi+ ls a 



wave- vector of HBdG ■ Evidently, Eqs. (45), (46) and (49) 
are, up to complications due to the sublattice, analogous 
to Eqs. (4) and (14) for the ID lattice. The sublat- 
tice degree of freedom really enters into the problem in 
the possibility of having inequivalent interfaces. This is 
made more formal by writing down the wave-matching 
conditions. For a type 1 boundary, they read 



SJ-l,m 
* S,l,m 



(50) 



with I = 0, L and for all m. For a type 2 boundary, 







( *SAm \ 


\IjB 

*e,£,m-l 




\I>2 













with I — 0, L and for all m. These four sets of equations 
define the transfer matrices, M^ } (0), M ( ^(L), M {2) (0) 

4 2) 



and M {2) (L). 



3. Decomposition of the transfer matrix and even odd effect 

We now express MA 1 '-') in terms of matrices in the unit- 
cell only. Using the transformation rule for a transla- 
tion of L unit-cells of length ax, we write 

jy[(*.i) — Diag[e~* fceiaT e ~ ik h La T e ik c La T e ik h La T ^ 

x M { 1\0) 



X Dl&g[e iksLaT e ~ ik s La T e ~ik s La T e ifcjLa T j 



xM«(0). 



(52) 



We introduce the following block notation for these ma- 
trices, 



a (3 
7 S 



•Mgj L) (0) 



3 L) • (53) 

v'fl(i) °R(L), 



To analyze the CAR processes, we have to compute 
t = (aT) _1 . We first factor out the free propagation 
by introducing 



a 



-ik E LaT 







-ikhLaT 



(54) 



(51) and using blockwise inversion formulae, we arrive at 



Jk s La T 



i2ksLa T 



1 







e -i2«He(fc s )£a T 



x [a 



(Oi 
n \ 



1 o \ w 

e i2y\e{k s )La T J 7fl 



(55) 



At this stage we ask for conditions (i) and (ii) to be 
true. As we said earlier, condition (ii), an even number 
of sites, is equivalent to the left-right symmetry i = j. 
Condition (i) implies in the case of graphene fRe^sa^) = 
or 2tt. Hence, for even lengths of the S region, 

with 

A=l 2 - [ai?]- l /3«[*W]-S« = 1 2 r«'r«\ (57) 

where and r& are two unitary reflection matrices, 
on the right and left interfaces respectively. Similarly 
to the ID case, we are left to prove that A has zero 



off-diagonal elements. To do so we try to find a simple 
relation between and . We take i = 1 for now, 

but a similar proof holds for i = 2. r^ 1 -*' and rjj are 
related by a translation of half a unit-cell pfots a sublattice 
exchange as is apparent in Fig. 10. We find 

r£ }/ = mag[-se ls ^ aT/2 e- is9 %-s'e- ls ' khaT/2 e is ' 0h }r^ ) 

xT>ia,g[-se lsk ' aT/2 e- ls(> %-s'e- ls ' khaT/2 e ls ' eh ], 

(58) 

for which a detailed proof is given in Appendix A. Due 
to particle-hole symmetry of the spectrum, the unitary 
matrix satisfies = a x r^a x . The off-diagonal 
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FIG. 10: r^ 3 ' and f£ )y are related by a translation of half a 

unit-cell. Due to the sublattice degree of freedom, f^ 13 ' and 

are not just related by a parity transformation but instead 
through phase factors (see Eq. (58)). 



elements of A are given by 



^ sg l e i(sk c a T /2-s'k h a T /2) e -i(sBe-s' 8 h ) 



(59) 



and 



Ahe _ Jl)he (l)hh ( -is' k h a 



t —2is8 h 



_\_ ss ' e i(sk c a T /2-s' k h a T /2) e -i(s0 c -s' 8 h ) \ ^gg-j 

Particle-hole symmetry implies k e — kh - see Eq. (47) - 
and in turn 8 e — 0^, as well as ss' = — 1. This leads to 
zero off-diagonal elements of A, and then of t. Hence, 
the CAR probability vanishes. 
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FIG. 11: (Color online) CAR probability in an anti-zigzag 
ribbon, for E F = V s = 0, A = 0.0031* and e = 0.003t, and 
W = 3 (number of chains). The inset shows the CAR prob- 
ability for a wider ribbon with W = 21, all other parameters 
kept equal. 



Qualitatively, electrons and holes in the lowest mode 
just propagate at the edges of the ribbon and therefore 
the situation is similar to the ID chain. The oscillations 
in the signal corresponding to the odd configurations is 
nevertheless unique to graphene with zigzag edges. 



IV. NON-IDEALITIES 
A. Armchair edges 



B. zigzag edges 

We now briefly turn to the case of graphene ribbons 
with zigzag edges. Whereas one can treat armchair rib- 
bons as quasi-lD wires, the inter-dependence of longitu- 
dinal and transverse momenta in zigzag ribbons makes 
such an approach difficult. We therefore computed the 
scattering matrix numerically, using a recursive Green's 
function algorithm. 36-40 We focus on the lowest mode, 
in order to simplify the analysis by avoiding inter-mode 
scattering. The main difference with respect to armchair 
graphene is a strong dependence on the parity of the 
width of the ribbon. For so-called zigzag edges (an even 
number of chains in the transverse direction) and zero 
Fermi energy, local (resp. crossed) Andreev reflection is 
suppressed in the lowest mode, due to opposite parities 
of the transverse wave-functions of incoming electrons 
and reflected (resp. transmitted) holes. 33 For anti-zigzag 
edges however, local Andreev reflection is allowed at zero 
Fermi energy, as incoming electrons and reflected holes 
have now the same parity. We found that for sharp NS 
interfaces, although incoming electrons and transmitted 
holes have still opposite parities, CAR is not necessarily 
suppressed since the pairing potential breaks the sublat- 
tice symmetry. In that case, an even-odd effect is indeed 
present, as shown in Fig. 11. 



The even-odd effect reported in the present paper re- 
lies on two stringent conditions, particle-hole symmetry 
of the spectrum, that is, E F = Vs = 0, and atom- 
ically sharp NS interfaces. There are several ways to 
relax these two conditions and we inspect them in turn, 
in the case of an armchair graphene nanoribbon. We 
use, for that purpose, two numerical methods. Either we 
compute elements of the scattering matrix by matching 
waves at the interface, numerically - this will be done 
for the simpler cases - or using a recursive Green's func- 
tions algorithm. 36-40 In all the examples we consider a 
metallic ribbon of width W — 92, and probe lengths of 
the S region up to a hundred ao, with ao the distance 
between carbon atoms, which amounts to a few nanome- 
ters. Energies are in units of t ~ 3eV, the nearest neigh- 
bor hopping. The superconducting gap is chosen to be 
A = Q.lt. For excitation energies below E g — 0.058£, 
only one mode, the metallic mode, is probed and con- 
tributes to the signal. This is the case in all plots in 
this section, where e — 0.03t. These are rather high val- 
ues, as compared to experiments. Indeed the induced 
gap is usually a few meV (~ 0.00K). 27-29 However the 
superconducting gates have typical lengths Lg of a few 
hundred nanometers. 27 ' 29 Therefore, with values realized 
in experiments, the dimensionless ratio ALs/hvF is of 
similar order as the same ratio computed with the values 
used in our calculation. 
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Note that, for a given length TV (measured in number 
of columns) of the S region, there are several possible 
configurations of the NS interfaces. We choose the 
following sweep: we fix the boundary in the unit-cell 
to be in configuration i = 1 and sweep the four 
configurations j = 1, 2, 3, 4 in the unit-cell L. 



Effects of doping 
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FIG. 12: (Color online) CAR probability for E F = O.lt, 
Vs — —Ef, A = O.lt and e — 0.03t, obtained by wave- 
function matching. Even and odd configurations are still 
clearly distinct. 

In this section we keep interfaces sharp. Let us look, for 
pedagogical purpose, at the following situation, Ep ^ 
and Vs = —Ep. This has two consequences. First, the 
real part of kg in the S region remains zero, hence we still 
expect a strong difference between even and odd config- 
urations, as can be inferred from Eq. (55). However, 
particle-hole symmetry is broken and the off-diagonal el- 
ements of the matrix A are not zero anymore. Then, the 
CAR probability does not vanish for even lengths of the 
S region, as can be seen in Fig. 12. 

If we now take Ep ^ 0, Vs = 0, that is, a constant 
doping throughout the sample, then the wave-vector ks 
acquires a non-zero real part, and the two curves start 
oscillating, as can be seen in Fig. 13. 

Next, we look at the case where there is only doping 
in the superconducting region, Ep = 0, Vs ^ 0. As can 
be seen in Fig. 14, the signal will, in general, show oscil- 
lations similar to the previous case of doping the whole 
sample. 

We nevertheless highlight a special situation, plotted 
in Fig. 15. The chosen value of Vs implies ^(ksap) — it. 
This has a very interesting consequence. A look at 
Eq. (55) allows us to infer that, for a given configura- 
tion of the interfaces - keep i and j fixed and vary L -, 
the signal does not oscillate with L, and in that sense, 
the situation is actually closer to that of Fig. 12. The 
fast oscillations are only between different configurations. 
However, the equivalence between configurations 1 and 3, 
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FIG. 13: (Color online) CAR probability for E F = O.lt, 
Vs = 0, A = O.lt and e = O.OSt, obtained by wave-function 
matching. 
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FIG. 14: (Color online) CAR probability for E F = 0, V s = 
O.li, A = O.lt, e = 0.03t (top panel) and E F = 0, V s = 
0.5t, A = O.lt, e = 0.03t (bottom panel). The red dashed 
line (resp. blue solid line) corresponds to even (resp. odd) 
configurations, and shows a long-wave length modulation due 
to doping. The signal oscillates fastly between the two curves 
(dotted line), a reminiscence of the even-odd effect. 



as well as between 2 and 4 is lifted. Indeed the CAR sig- 
nal clearly interpolates between four different curves in 
Fig. 15, contrary to, for instance, the case of Fig. 14. 



The transmission for the configuration (1,3) can be de- 
duced from Eq. (55) by taking i = j = 1 but replacing 
L by L + 1/2 (remember that 1 and 3 are related by a 
translation of half the unit-cell). Now if Dlz(ksaT) was 
strictly 7r, the matrix in Eq. (55), for the configuration 
(1,3) would be a z and not the identity as for the configu- 
ration (1,1). This was exactly what was happening in the 
ID case, as one can see in Eqs. (25) and (26). Of course, 
the CAR does not vanish exactly because particle-hole 
symmetry is obviously broken, however the signal shows 
two intertwined very strong oscillations, as visible in the 
inset of Fig. 15. 
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FIG. 15: (Color online) CAR probability for E F = 0, 
V S = 1.41t, A = O.lt, e = O.OSt. For this value of Vs, 
9ie(ksaT) — 7!"- The CAR signal (dotted line) interpolates 
between four curves: Long- wave length oscillations are sup- 
pressed and configurations (1,1) - circles - and (1,3) - rhombi 
- as well as (1,2) - squares - and (1,4) - triangles - cease to 
be equivalent. The inset shows the same signal, highlighting 
the two intertwined fast oscillations (see text). 



2. Effects of smooth NS interfaces 

We now investigate the effect of smooth NS inter- 
faces. Each plot in this section corresponds to a situ- 
ation where the pairing potential goes smoothly from its 
bulk value A to zero in the normal leads over a length 
da (resp. cIl) on the right interface (resp. left inter- 
face), as sketched in Fig. 16. The smoothing function 
is defined as 0.5 (cos((y — N)/djf) + 1) on the right in- 
terface and 0.5 (cos(y/di) + 1) on the left one. dji, and 
dz, are of the order of the coherence length in the su- 
perconductor, £ = hvp/A. The signal is plotted against 
N, the length of the superconducting region over which 
the gap is constant (see the inset of Fig. 16). We have 
made a numerical calculation of the scattering matrix, 
using a recursive Green's function algorithm 36-40 . We 
highlight a few possible effects. At zero Fermi energy the 
CAR probability is strongly suppressed, as can be seen in 
Figs. 16, 17 and 18. Indeed, as one makes the NS inter- 
faces smoother, the system should be better described by 
the Dirac-Bogoliubov-de Gennes equations with a pairing 
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potential smooth at the level of the sublattice. In that 
case, CAR would be completely suppressed at zero Ep, 
due to the orthogonality of the spinor wave-functions of 
incoming electrons and reflected holes. However, we find 
that the even-odd effect is preserved in the case of a sym- 
metric smoothing, that is dji — dp and d^ is commen- 
surate with the length of the unit-cell, as can be seen in 
Fig. 16. We obtain that the zeros of the CAR probability 
are lifted by breaking the left-right symmetry of the NSN 
junction, which can be done either by taking da ^= d^, 
as in Fig. 17, or by choosing a smoothing length that is 
incommensurate with the unit-cell length, as in Fig. 18. 
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FIG. 16: (Color online) CAR probability for E F = 0, V s = 0, 
A = O.lt, e = 0.03t. The superconducting potential is 
smoothed over a distance dn = g(l = 9an - with an the dis- 
tance between carbon atoms - corresponding to three unit- 
cells. The even-odd effect is preserved. 
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FIG. 17: (Color online) CAR probability for E F = 0, V s = 0, 
A = O.lt, e = 0.03t, d L = 9a , and d R = 18a . The left- 
right symmetry of the junction is broken. The signal (dotted 
black line) oscillates between four configurations showed with 
different symbols (squares, circles, rhombi and triangles). 

These two possibilities share a common feature: As 
one increases the length of the S region, the signal 
interpolates between four distinct curves, signaling four 
inequivalent configurations. 

Next, we combine smoothing and doping by consid- 
ering the case of finite doping in the S region, while 
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FIG. 18: (Color online) CAR probability for E F = 0, V s = 0, 
A = O.lt, e = 0.03i, d L = 10oo, and d R = Wa . The left-right 
symmetry of the junction is also broken, as the smoothing 
length is incommensurate with the unit-cell length. 



preserving the left-right symmetry. We picked two dif- 
ferent values of doping, Vs — 0.1 and Vs — 0.5. The 
corresponding Figs. 19 and 20 are to be contrasted with 
Fig. 14. Note that the Fermi energy goes smoothly from 
in the S region to in the normal leads on the same 
distance than the pairing gap. Finite doping seems to 
smooth out every trace of the even-odd effect as there 
is hardly any difference between possible configurations 
of the interfaces. In Fig. 19, the four curves (obtained 
by increasing the length of the S region by one unit-cell 
starting from the four a priori different configurations) 
are still distinct but their interpolation is so smooth that 
the sublattice hardly matters. This is even clearer for the 
larger value of doping in Fig. 20, where one cannot dis- 
tinguish between any sublattice configurations, and the 
system is better described by the Dirac equation. One 
can check that by decreasing the smoothing length, the 
oscillations due to the lattice reappear, as shown in the 
inset of Figs. 19 and 20. 
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FIG. 19: (Color online) CAR probability for E F = 0, V s = 
O.lt, A = O.lt, e = 0.03t, d L = 9a , and d R = 9o - The 
four configurations that were visible in Figs. 17 and 18 now 
lie practically on top of each other. In the inset we have 
taken d R — dL — 3ao. Fast oscillations between different 
configurations are more visible. 



FIG. 20: (Color online) CAR probability for E F = 0, V s = 
0.5*, A = O.lt, e = 0.03i, d L = 9a , and d R = 9a . We still 
use different symbols for the four configurations of Figs. 17 
and 18, although the signal interpolates smoothly between 
them (see text). In the inset we have taken d R = d_L = 3ao. 
Fast oscillations between different configurations are more vis- 
ible. 
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FIG. 21: (Color online) CAR probability for E F = 0, V s = 
lAlt, A = O.lt, e = 0.03t, d L = 9a , and d R = 9a . For 
this value of Vs, 9\e(ksa,T) — it- The four configurations, 
showed with different symbols, are again clearly visible. The 
inset is the same signal, highlighting the two intertwined fast 
oscillations (see text). 

However, the special case highlighted in Fig. 15 is not 
altered by smoothing - although the amplitude of the 
signal is damped - as is shown in Fig. 21. This is a case 
where oscillations strictly due to doping are suppressed 
and fast oscillations due to the positions of the interfaces 
with respect to the sublattice re-emerge. We still identify 
two intertwined very fast oscillations, reminiscent of the 
even-odd effect. 



B. zigzag edges 

We now briefly turn to ribbons with zigzag edges 
and show how the even-odd effect is altered by doping 
and smoothing. The effect of doping is similar to the 
armchair case: The CAR probability exhibits long wave- 
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length oscillations, yet the even and odd configurations 
remain clearly visible, as shown in Fig. 22. 
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FIG. 22: (Color online) CAR probability in an anti-zigzag 
ribbon, for E F = 0, Vs = O.OOli, A = 0.0031i and e = 0.003i, 
W = 3 (number of chains) . The inset shows the CAR prob- 
ability for a wider ribbon with W = 21, all other parameters 
kept equal. 

The effect of smoothing however is quite different and 
depends strongly on the width of the ribbon, as shown 
in Fig. 23. Similar to the armchair case, we find that 
the CAR probability is strongly suppressed at zero Fermi 
energy. The even-odd effect is preserved for narrower 
ribbons while it is lifted as the width increases. 
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FIG. 23: (Color online) CAR probability in an anti-zigzag 
ribbon, for E F = = Vs = 0, A = 0.0031t and e = 0.003i, 
W = 3 (number of chains), and dn = d,L = 20\/3fflo (20 
unit-cells). The inset shows the CAR probability for a wider 
ribbon with W = 21, all other parameters kept equal. 



V. CONCLUSION 

We have investigated the variations of the CAR prob- 
ability as a function of the length of the superconduct- 
ing region, and showed that, in a variety of situations, 
large amplitude oscillations on the scale of the lattice 
are present. We have discovered conditions under which 



the CAR probability vanishes altogether only for even 
lengths of the S region, and explained this effect on the 
basis of symmetries of the transfer matrix. These suffi- 
cient conditions are not strictly tied to a specific lattice 
and we found that they are satisfied at least in both the 
simple ID lattice and graphene nanoribbons. 

Although perturbations tend to lift the perfect zeroes 
- for instance, by breaking the left-right symmetry 
or the particle-hole symmetry - the CAR amplitude 
often retains fast oscillations, a remnant of the even- 
odd effect. It is worth emphasizing that, although 
our analytical proof of the effect relies on sharp NS 
interfaces, smoothing of the interfaces, as long as the 
left-right symmetry is not broken, preserves the effect, 
as we have shown numerically for graphene. However, 
the amplitude of the signal is very much reduced, as 
one approaches the Dirac regime, where selection rules 
would constrain the CAR probability to be exactly 
zero in the absence of doping. When doping in the S 
region is added to the smooth boundaries, in armchair 
nanoribbons the damping is reduced and fast oscillations 
tend to disappear, as the smoothing distance approaches 
and exceeds the coherence length. However, we exhibit a 
point of very high doping in the S region where the fast 
oscillations are barely affected by smoothing. Although 
there is no strict even-odd effect in that case - the 
CAR probability does not vanish exactly - the large 
amplitude of the oscillations is readily understood on 
the basis of the transfer matrix. It makes yet another 
connection with the ID lattice, and hints at the possibil- 
ity of finding similar effects in seemingly different lattices. 
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Appendix A: Transfer matrices and wave-matching 

1. General properties of transfer matrices 

The transfer matrix and the scattering matrix of an 
NSN junction are 4x4 matrices of the form 



r t' 



M 



t r' / ' V 7 5 1 ' 



OL (3 



(Al) 



where r, t, t', r' are the usual reflection and transmission 
matrices, from both sides of the barrier. These 2x2 
matrices are related to the blocks of M through 41 



r = -5- 1 1 , 
t = (at)" 1 , 



r'=(3S-\ 
t' = d-\ 



(A2) 
(A3) 
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In the case of time-reversal symmetry, then S — a* and 
7 = (3* . In Eq. (28), we have introduced the reflection 
matrices r# and r' L at the NS interfaces of Fig. 6. From 
Eq. (A2) we have tr = —d^fn and r' L = PlS^ • Using 
inversion formulas and Mi(0) = [M/f (0)] _1 it follows 
that r' L = —a^fln. This proof also holds for armchair 
graphene nanoribbons. 



2. Important relations in the case of the ID lattice 

We write down explicitly the wave-matching equations 
for the ID lattice and prove important relations between 



Ol,e + h,e 

a\,h + hfi 

a he e ik ° a + b he e- ik ° a 
ai, h e~ lkha + b hh e lk » a 



where F+ was defined in Eq. (1 



&2,e + a 2 , e 
h.h + a2,h 

b 2<e e ik ° a + a 2 . e e~ ik ' a 



as is clear from Fig. 5. Inverting these two systems we 
would get 

(X)= m «o>(b;)- < a6 > 

and 

respectively. However one can check by inspection that, 
for given and Ag, B 2 = Ai and A 2 = Bi are still 
solutions of the system (A5). Therefore we have, 

(b;)=^<»>(a s ,). < as > 

and, as a consequence, M. L (0) — [M fl (0)] _1 . 



transfer matrices. 



Let us first prove Eq. (22), M L (0) = [M^O)] -1 . 
The transfer matrix Mft(O) is defined through the wave- 
matching equations in Eq. (16), which explicitly read, 
for e > 0, 



I 

We now prove Eq. (29). To do so one should write 
the system of wave-matching equations for the transfer 
matrix Mt(0) defined in Fig. 6. Comparing it to the 
system defining Mfl(0), one finds the simple relation 

^(°) = (i° 2 o) [ M ^" (i 2 o)' (A9) 

from which follows r' L = r# and then Eq. (29). 

3. Important relations in the case of armchair 
graphene nanoribbons 

The proof of M^(0) = [M^(O)]- 1 is straightforward 
and utterly similar to the proof for the linear chain. In 
the present section, we focus on the proof of Eq. (58), 
which comes from a relation very similar to Eq. (A9). 
Again we prove the relation by simple inspection. Let 



=b s ,iF+(k s ) + b s ,2F + (k s )* + a s ,iF+(k s ) + a s ^F + (k s )* , 
=bs,i + b.s,2 + as,i + asa , 

=b sa F + (k s )e iksa + b s , 2 F + (k s y e - lk > + a s ,iF + (k s )e- lksa + as^+fe) V fc * a , 
=b s ,ie lksa + b S ae~ tk * sa + a s ,ie-* ksa + a s , 2 e lk s a , (A4) 



I 

2). Similarly, Mi(0) is defined through the equations 



=b SA F+(k s ) + bs, 2 F + (ksT + as,iF+(k s ) + a s , 2 F + {k s )* 
=bs,x + b s ,2 + as,i + as,2 

=b s .iF + (k s y k * a + b s . 2 F + (k s )*e- lk s a + a SA F + (k s )e- iksa + a s , 2 F + (k s )* e lk > 
=b s ,ie iksa + b s , 2 e~ lk s a + a s ,ie-' lksa + a s , 2 e ik ' sa (A5) 
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us write the wave-matching equations for both M. ( ^'(0) the scattering states defined in Appendix B for graphene. 
and M^(0), as defined in Fig. 10. To that end we use First for M£(Q), in the case E F = V S = 0, 



a 1 , e {-ae- u0 -) + h^-se™ 6 *) =b s ,i(-e- ies )F + b s , 2 e i6s F* + a S|1 (-e* s )F + as^'F* 

ai ^ e -isKa t /2 + biee isk e a T /2 =b Sil Fe- iksaT/2 + b Sy2 F* e ~ lksaT I 2 + a s ,iFe lksaT/2 + a s , 2 F* e iksaT/2 

a hh {-s'e is ' e ») + b 1>h (-s'e- is ' e ») =&s,i(-e^ s ) + b s ^ 9s + a S)1 (-e i9s ) + ag^ 3 

aihe ™'k h a t /2 + biee -ts'k h a T /2 =&gie -ifc S ax/2 + ^g-ifcso-r/2 + ^^ikgar/i + as ^k s a T /2 (A10) 



where we defined F = F±+ to simplify notations. For 



M^(0), still in the case E F = V s = 0, we have 



b 2 ,e{-se- ls6 ')e lsk " at/2 + a 2 ^{-se ls9 ")e- isKat / 2 =a SA {-e~ Ws )Fe iksaT I 2 + a Si2 e Ws F* e iksaT/2 

+ b s ,i(-e if)s )Fe- iksaT/2 + b sa e~ i9s F*e~ lksaT/2 
b 2 , e + a 2 ,e =as.\F + a s , 2 F* + bs,\F + b s , 2 F* 

b 2 , h {-s'e is ' eh )e- ls ' khat/2 + a 2 M(-s , e^ s ' 9h )e +ls ' khaT/2 =a s ,i{-e- l9s )e lksaT /2 + a s , 2 e i9s ' e iksaT/2 + b s ^{-e l9s )e~ lksaT /2 

+ b s , 2 e- l9s e- lksaT/2 

b 2 ,h + a 2 ,h =asA + a s , 2 + b s ,i +b s , 2 (AH) 



It is just a matter of algebra to verify that the transfor- 
mation, 



Ai 



1 2 
1 2 



and 



B s 
As 



Diag 



1 2 
1 2 



-se- ls9 " , -s'e is 9h , -se ts9 " , -s'e~ is 9h 



B 2 

A 2 



Diag [- 



(A12) 



(A13) 



maps the system (A10) to the system (All), from which 
we deduce 



Mi 1) (0)=P eh l i ; ^1 |M^(0) 



1. 



1 2 
1 2 



V S , 
(A14) 



where we have defined 
Veh = Diag 



se- is9 %-s'e ls ' 9 \-se ls9 '-s'e- ls ' 9h 



(A15) 



and 

V s = Diag [-e l9s , e~ i8s , -e~ l9s , e %9s ] , (A16) 
with 9s = 0(ks)- Equation (58) follows directly. 

Appendix B: Scattering states in graphene 

In this appendix, we detail the solutions of the 
Bogoliubov-de Gennes equations in both normal and su- 
perconducting regions. In the normal regions, right (+) 
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and left (— ) going waves are defined for electron-like ex- 
citations as 



with 



,(*) = (<P 



± „±isk c yt, m A 



± ±isk c ye t „ 



—se T 
1 



( B1 ) 



(B2) 



s = sgn(e + Ep) and 8 e = 9(k e ). We have dropped 
the sine function, since it simplifies in the wave-matching 
conditions. For hole-like excitations 



(B3) 



with 



4A = f- s ' e±ls ' 6h 

4b) \ i 



(B4) 



sgn(— e + Ep), and 6^ = 9(kf l ). The wave-vectors k e Diag[e 



and kh are the positive roots of e = —Ep + E s (k,p) and 
e = Ep — E s r(k,p), respectively. Particle-hole symmetry 
of the spectrum, Ep = 0, implies 



k e — Ai/j. 



(B5) 



In the superconducting region, where A =/= 0, 
the BdG Hamiltonian has four eigenva lues labeled 
ei± = ± y /A 1 (k) 2 + A 2 and e 2± = ±y / A 2 (k) 2 + A 2 
with Ai(fc) = E + (k,p) - E F - V s and A 2 (fc) = 
E + (k,p) + Ep + Vs- The corresponding eigenvectors 

are Xl± (fc) = {-e~ ie ^ F l± {k), F 1± (k), -e~^ k \ if 
and X 2±(k) = (e- i6 ^ F 2± (k),F 2± ( k),e- i6 ^,l) T with 
Fi±(k) = (Ai(fc) ± v/AT(^F+A2)/A and F 2± (fc) = 
(-A 2 (fc) ± v/A 2 (A:) 2 + A 2 )/A. As for the linear chain, 
there are four available wave-vectors, at a given excita- 
tion energy. However, we must distinguish between the 



cases, Ep + Vs ^ and Ep + Vs — 0. Indeed, in the 
former case, all modes stem from one eigenvalue, e.g. £1+ 
for e > and Ep + Vs > 0. For e < A, the wave-vectors 
are given by 



* . ± A Arccos (WjgEgd-tj) . 

a T \ 2e p J 

(B6) 

We write these wave- vectors as fcg, — fcg (right-decaying 
modes) and — ks, k* s (left-decaying modes) where ks is 
the solution with positive real and imaginary parts. A 
general solution of the BdG equations in the S region, for 
< £ < A and Ep + Vs > is a superposition of right 
and left decaying modes of the form 

*s,e,rn = b s ,i®(ks)xi+(ks) + bs,2®{-k s )xi+(-k*s) 

+ a s ,i®(-k s )xi+(-ks) + a s ,2®(k s )Xi+(k s ), 

(B7) 

where we have introduced the diagonal matrix D(k) — 



iky, 



iky e 



,ikyt 



ikyr 



In the case 



Ep + Vs = 0, ei+ and e 2+ (rcsp. ei_ and e 2 _) are degen- 
erate. There are still four modes but only two different 
wave-vectors: ks and — ks, with 3m(ks) > 0. Then, 
a general solution in the S region, for < e < A and 
Ep + Vs = is of the form 

^S,i,m = b s ,i'Zi{ks)xi+{ks) + bs,2®{k s )x2+{k s ) 

+ a s ,i®(-k s )xi+(-ks) + a s .2®(-ks)x2+(-k s ). 

(B8) 

We emphasize the important constraint from particle- 
hole symmetry of the spectrum that is 



mc(ks) 



if e p < 0, 

2n/aT if e p > 0. 



(B9) 
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